#=============================================================================#
#
# PROJECT:        Does Islamist Terrorism Still Affect Political Attitudes?
# AUTHORS:        ** anonymized for review **
# CONTACT:        ** anonymized for review **
# LAST MODIFIED:  September 15, 2025
# 
#=============================================================================#
#
# This R file contains the code used to replicate all results reported in the 
###### main paper.
#
##### output:
# Table 1
# Table 2
# Figure 1
# 
#=============================================================================#

rm(list=ls())
getwd()

# Install and load all necessary packages with the ipak function:
    # i.e., check to see if packages are installed. 
    # Install them if they are not, 
    # then load them into the R session.

ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

packages <- c("readxl", "metafor", "metaSEM", "tidyverse",
              "ggplot2", "modelsummary", "stargazer") 
ipak(packages)



# Data ------------------------------------------------------------------------

df <- readRDS("data/df.rds") 
# delete 3 studies on effects of terrorism conducted outside Western democracies 
# extra criterion 3
df <- df %>% 
  filter((ExactAttack_2 != "2002 Bali bombings" &  # attacks on western countries 
            ExactAttack_2 != "2008 Mumbai attacks" &  # attacks on western countries 
            ExactAttack_2 != "Gaza - Israel conflict") %>%
           replace_na(TRUE))
# summary
summary(df)


# Main Analyses ----------------------------------------------------------------

## Average effect size --------------------------------------------------------

## Get estimates for...
# ... Overall sample
m0 <- rma.mv(d, var_d, 
             random = ~ 1 | ID_R/ID_ES_Unique, 
             data = df)
m0


## Average effect size over time ----------------------------------------------

df %>% count(is.na(StudyYear)) %>% mutate( prop = n / n_distinct(df))
df %>% filter(!is.na(StudyYear)) %>% distinct(ID_R) %>% nrow()
df <- df %>%
  mutate(StudyYear = StudyYear - 2001)
df %>% 
  summarize(
    min = min(StudyYear, na.rm = TRUE),
    max = max(StudyYear, na.rm = TRUE)
  ) # study year: range


## Get estimates for...
# ... Model 1
m1 <- rma.mv(d, var_d, 
             mods = ~ StudyYear, # StudyYear as moderator
             random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
             data = df)
m1

# ... Model 2
m2 <- rma.mv(d, var_d, 
             mods = ~ StudyYear + # StudyYear as moderator +
               CasualtiesLess10 + Casualties10to100 + Casualties100plus +
               Exp + NatExp +
               ProbSample +
               Outgroup + Rally +
               US,
             random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
             data = df)
m2


# Explore the correlation between Year of Study and being US-based
df$US <- as.numeric(df$US) - 1
mean(df$US[df$StudyYear + 2001 < 2005], na.rm = TRUE)
mean(df$US[df$StudyYear + 2001 >= 2005], na.rm = TRUE)
cor.test(df$StudyYear, df$US)

summary(df)
library('rcompanion')
cramerV(df$US, df$NineEleven, bias.correct = FALSE)
cor.test(as.numeric(df$US), as.numeric(df$NineEleven))


## Table (Table 1)

# Step 1: Create a named list of models
models <- list(
  "Model 1" = m1,
  "Model 2" = m2
)

# Create the table
cm <- c("StudyYear"     = "Year of study",
        "CasualtiesLess101"  = "Fatalities: $<$10",
        "Casualties10to1001" = "Fatalities: 10–100",
        "Casualties100plus1" = "Fatalities: $>$100",
        "Exp1"          = "Randomized experiment",
        "NatExp1"       = "Natural experiment",
        "ProbSample1"   = "Probability sample",
        "Outgroup"      = "Outgroup hostility",
        "Rally"         = "Rally tendencies",
        "US1"           = "US dummy",
        "intercept"     = "Constant")
modelsummary(models, 
             coef_map = cm,
             stars = c("*" = .05, "**" = .01, "***" = .001),
             fmt = 3,
             gof_omit = 'IC|Log|Adj|AIC|BIC|RMSE',
             escape = FALSE,
             output = "latex_tabular")



## Effect of past attacks on effect sizes -------------------------------------

## Get estimates for...
# ... Model 2: Attacks in the same country during the 5 years prior to the study.
m2_past_1 <- rma.mv(d, var_d, 
                    mods = ~ PastAttacks_Country_5yr + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      Outgroup + Rally +
                      US,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_1

# ... Model 2: Attacks in the West during the 5 years prior to the study.
m2_past_2 <- rma.mv(d, var_d, 
                    mods = ~ PastAttacks_West_5yr + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      Outgroup + Rally +
                      US,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_2

# ... Model 2: Attacks in the same country, no time restriction.
m2_past_3 <- rma.mv(d, var_d, 
                    mods = ~ PastAttacks_Country_All + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      Outgroup + Rally +
                      US,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_3

# ... Model 2: Attacks in the West, no time restriction.
m2_past_4 <- rma.mv(d, var_d, 
                    mods = ~ PastAttacks_West_All + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      Outgroup + Rally +
                      US,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_4


## Table (Table 1)

# Step 1: Create a named list of models
models <- list(
  "Model 1" = m2_past_1,
  "Model 2" = m2_past_2,
  "Model 3" = m2_past_3,
  "Model 4" = m2_past_4
)

# Create the table
cm <- c("PastAttacks_Country_5yr" = "Number of previous Islamist attacks (same country, 5 years)",
        "PastAttacks_West_5yr" = "Number of previous Islamist attacks (all Western countries, 5 years)",
        "PastAttacks_Country_All" = "Number of previous Islamist attacks (same country, all years)",
        "PastAttacks_West_All" = "Number of previous Islamist attacks (all Western countries, all years)",
        "intercept"     = "Constant")
modelsummary(models,
             coef_map = cm,
             stars = c("*" = .05, "**" = .01, "***" = .001),
             fmt = 3,
             gof_omit = 'IC|Log|Adj|AIC|BIC|RMSE',
             escape = FALSE,
             output = "latex_tabular")


## Average effect size by attack ----------------------------------------------
df_attacks <- df %>% 
  filter(Attack == 1, # reference to/study real attacks
         ExactAttack_2 != "2014 Jewish Museum of Belgium shooting", 
         # only 1 observation: no Meta-Analysis possible
         ExactAttack_2 != "2015 San Bernardino",
         # only 1 observation: no Meta-Analysis possible
         ExactAttack_2 != "Hypothetical attack", # no hypothetical attacks
         ExactAttack_2 != "Other attack or combination", # no combinations
         ExactAttack_2 != "Unspecified attack", # no unspecified attacks
         ExactAttack_2 != "2002 Bali bombings", # attacks in western countries 
         ExactAttack_2 != "2008 Mumbai attacks", # attacks in western countries 
         ExactAttack_2 != "Gaza - Israel conflict") # attacks on western countries 

df_attacks %>% distinct(ExactAttack_2) %>% nrow()
unique(df_attacks$ExactAttack_2)
df_attacks %>% distinct(ID_ES_Unique) %>% nrow()
df_attacks %>% distinct(ID_R) %>% nrow()


## Get estimates for...
# ... Model 0, with accounting for attack-clustering
df_attacks$ExactAttack <- as.factor(df_attacks$ExactAttack_2)
m0_attacks <- rma.mv(d, var_d, 
                     random = ~ 1 | ExactAttack/ID_R/ID_ES_Unique, # four-level model
                    data = df_attacks)
m0_attacks

### models
n.attacks <- length(unique(df_attacks$ExactAttack_2))
my.cool.loop.results <- data.frame( attack = rep(NA, n.attacks),
                                    n.k = rep(NA, n.attacks),
                                    n.j = rep(NA, n.attacks),
                                    coef = rep(NA, n.attacks),
                                    lbound = rep(NA, n.attacks),
                                    ubound = rep(NA, n.attacks) )
for(i in 1:n.attacks) {
  my.attack <- unique(df_attacks$ExactAttack_2)[i]
  meta_attacks <- rma.mv(d, var_d,
                         random = ~ 1 | ID_R/ID_ES_Unique,
                         data = df_attacks[df_attacks$ExactAttack_2==my.attack,])
  # start filling in dataframe
  my.cool.loop.results$attack[i] <- unique(df_attacks$ExactAttack_2)[i]
  my.cool.loop.results$n.k[i]   <- summary(meta_attacks)$s.nlevels[1]
  my.cool.loop.results$n.j[i]   <- summary(meta_attacks)$s.nlevels[2]
  my.cool.loop.results$coef[i]  <- summary(meta_attacks)$b[1]
  my.cool.loop.results$lbound[i]   <- summary(meta_attacks)$ci.lb
  my.cool.loop.results$ubound[i]   <- summary(meta_attacks)$ci.ub
}

# results ordered by year/attack (for plot)
my.cool.loop.results <- my.cool.loop.results[order(my.cool.loop.results[,1],decreasing=FALSE),]
my.cool.loop.results 
stargazer(my.cool.loop.results, summary=FALSE)

# results ordered by size (in text)
my.cool.loop.results[order(my.cool.loop.results[,4],decreasing=FALSE),] 


### figure
graphics.off()

# get estimates and confidence intervals
d <- matrix(NA,length(my.cool.loop.results$attack),1)
ci <- matrix(NA,length(my.cool.loop.results$attack),2)

for (i in 1:length(my.cool.loop.results$attack))
{
  d[i] <- my.cool.loop.results[i,]$coef
  ci[i,] <- cbind(my.cool.loop.results[i,]$lbound,
                  my.cool.loop.results[i,]$ubound)
}

# safe names of attacks + n's
attacks <- c(expression("9/11"*" (N"[j]*"=33, N"["i"]*"=178)"),
             expression("Madrid 2004"*" (N"[j]*"=9, N"["i"]*"=43)"),
             expression("van Gogh 2004"*" (N"[j]*"=5, N"["i"]*"=32)"),
             expression("London 2005"*" (N"[j]*"=9, N"["i"]*"=56)"),
             expression("Fort Dix 2007"*" (N"[j]*"=1, N"["i"]*"=3)"),
             expression("Boston 2013"*" (N"[j]*"=3, N"["i"]*"=6)"),
             expression("Ottawa 2014"*" (N"[j]*"=1, N"["i"]*"=13)"),
             expression("Paris 2015a"*" (N"[j]*"=10, N"["i"]*"=56)"),
             expression("Copenhagen 2015"*" (N"[j]*"=1, N"["i"]*"=5)"),
             expression("Paris 2015b"*" (N"[j]*"=19, N"["i"]*"=80)"),
             expression("Berlin 2016"*" (N"[j]*"=5, N"["i"]*"=39)"),
             expression("Brussels 2016"*" (N"[j]*"=4, N"["i"]*"=50)"),
             expression("Manchester 2017"*" (N"[j]*"=1, N"["i"]*"=2)"),
             expression("Stockholm 2017"*" (N"[j]*"=1, N"["i"]*"=12)"),
             expression("Strasbourg 2018"*" (N"[j]*"=1, N"["i"]*"=5)"))

# plot results
par(mar = c(10, 4.5, 1.5, 0.1))  

y.lab <- seq(from=1, to=length(my.cool.loop.results$coef), by=1)
x.lab <- seq(from=1, to=length(my.cool.loop.results$attack), by=1)

plot(x.lab, y.lab,
     type="n",  bty = "n", las = 1,
     ylab = "", xlab = "", 
     xaxt="n", 
     ylim = c(-0.2,0.6),
     cex.lab=1.5, cex.axis=1.3, cex.main=1.5)
axis(1, at=x.lab, labels=F)
text(x=x.lab, y=par()$usr[3]-0.03*(par()$usr[4]-par()$usr[3]),
     labels=attacks, srt=45, adj=1, xpd=TRUE, cex = 1.1) +
mtext(expression(paste("Average effect sizes (Cohen's ", italic("d"), ")")),
        cex = 1.2, side = 2, line = 3, outer = F)
points(x.lab, d, pch = 19, col = "black")
segments(x.lab, ci[,1], x.lab, ci[,2], lty = 1, lwd = 2, col = "black")
abline(h=0, col=1, lwd=2, lty = 1)

abline(h=0.1612, col=8, lwd=2, lty = 2)
text(12, 0.185, col=8, "Average effect of attacks") 
#arrows(6.5, 0.21, 8, 0.50, col=8, lwd=2, lty = 2)

# safe plot
meta.attacks <- recordPlot()
pdf(file="output/figures/meta_attacks.pdf",
    width = 10,
    height = 7)
meta.attacks
dev.off()
